# author: rehmertjochen@gmail.com
# code to replicate:
# Tables:  1, 5, 7 - 15
# Figures: 4, 5, 12

rm(list=ls())
# directory
setwd("C:/Users/rehmert.PWIREHMERT-NB/Dropbox/GreenSelection/replication")

# data
library(lmtest)
library(multiwayvcov)
library(sandwich)
library(stargazer);library(clubSandwich)
dat <- readRDS("incumbents.RDS")


# function for obtaining average marginal effects
ame.p <- function(input, dv = c("f_denied","rank_diff","rank_final_prob", "ranks_competed","numcand_final", "margin_final"),
                  iv = c("theta_mean_land_deviation","pq_both_coverage_land","pq_both_num_land"), 
                  seed = 1104, nsim = 1000){
  require(sandwich);require(lmtest);require(MASS)
  dat.tmp <- input[["model"]]
  colnames(dat.tmp)[which(colnames(dat.tmp) %in% "as.factor(incumbent)")] <- "incumbent"
  if(dv %in% c("rank_diff","margin_final") ){form1 <- as.formula(paste(dv,"~", paste(names(dat.tmp[-1]), collapse ="+") ))}else{form1 <- input[["formula"]]}
  #if(dv %in% c("rank_final_prob") ){form1 <- as.formula(paste(dv,"~", paste(names(dat.tmp[-1])[-which(names(dat.tmp[-1]) == "smd_nomination")], collapse ="+") ))}else{form1 <- input[["formula"]]}
  
  # run models
  if(dv == "f_denied"){mod1 <- glm(form1, data = dat.tmp, family = binomial("logit"))}
  #if(dv == "rank_final_prob"){mod1 <- glm(form1, data = dat.tmp, family = quasibinomial("logit"))}
  if(dv %in% c("rank_diff","margin_final")){mod1 <- lm(form1, data = dat.tmp)}
  if(dv == "ranks_competed"|dv == "numcand_final"){mod1 <- glm(form1, data = dat.tmp, family = poisson("log"))}
  
  require(fastDummies)
  mod.dat.0 <- mod.dat.1 <- data.frame(cbind(dummy_cols(dat.tmp$incumbent)[-1], mod1[["model"]][-c(1:2)]))
  mod.dat.1[[iv]] <- dat.tmp[[iv]] + mean(dat.tmp[[iv]], na.rm = TRUE)
  
  # simulations mod1 cov(mod1)
  set.seed(seed)
  
  b.sim1 <- mvrnorm(nsim, coef(mod1), vcovHC(mod1, type="HC"))
  
  lin.preds.0 <- as.matrix(mod.dat.0) %*% t(b.sim1)  
  lin.preds.1 <- as.matrix(mod.dat.1) %*% t(b.sim1)  
  
  if(dv %in% c("f_denied")){   # , "rank_final_prob"
    preds.0 <- exp(lin.preds.0)/(1+exp(lin.preds.0))
    preds.1 <- exp(lin.preds.1)/(1+exp(lin.preds.1))
  }
  if(dv %in% c("ranks_competed","numcand_final")){
    preds.0 <- exp(lin.preds.0)
    preds.1 <- exp(lin.preds.1)
  }
  if(dv %in% c("rank_diff","margin_final")){
    preds.0 <- lin.preds.0
    preds.1 <- lin.preds.1
  }
  
  preds.0.mean <- apply(preds.0, 2, mean, na.rm = TRUE)
  preds.1.mean <- apply(preds.1, 2, mean, na.rm = TRUE)
  
  ame.95 <- quantile(preds.1.mean - preds.0.mean, probs = c(0.5, 0.025, 0.975), na.rm = TRUE)
  
  ame.90 <- quantile(preds.1.mean - preds.0.mean, probs = c(0.5, 0.05, 0.95), na.rm = TRUE)
  
  return(list(iv = iv,
              dv = dv,
              ame90 = ame.90,
              ame95 = ame.95))
  
}



#################################################
# 1) number of ranks competed
#################################################
dat$ranks_competed[!is.na(dat$ranks_competed)] <- dat$ranks_competed[!is.na(dat$ranks_competed)] -1

library(AER)
##### -- All -- #####
### --- PQs - Share --- ###
summary(a2.c <- glm(ranks_competed ~ - 1+ as.factor(incumbent) + rank_before_prob + 
                      office + elected+ theta_mean_land_deviation + ln_pq_total+ pq_both_coverage_land ,  data = dat, family = poisson("log")))
a2.c.se.rob <- coeftest(a2.c, vcov = vcovHC(a2.c, "HC"))[,2];coeftest(a2.c, vcov = vcovHC(a2.c, "HC"))
a2.c.se.cl <- coef_test(a2.c, vcov = "CR2", 
                        cluster = a2.c[["model"]]$`as.factor(incumbent)`, test = "Satterthwaite")[,"SE"]
rank.cov.a <- ame.p(a2.c, dv = "ranks_competed", iv = "pq_both_coverage_land",seed = 1104, nsim = 1000)

### --- PQs - Number --- ###
summary(a2.d <- glm(ranks_competed ~ - 1+ as.factor(incumbent) + rank_before_prob + 
                      office + elected+ theta_mean_land_deviation + ln_pq_total+ ln_pq_land ,  data = dat, family = poisson("log")))
a2.d.se.rob <- coeftest(a2.d, vcov = vcovHC(a2.d, "HC"))[,2];coeftest(a2.d, vcov = vcovHC(a2.d, "HC"))
a2.d.se.cl <- coef_test(a2.d, vcov = "CR2", 
                        cluster = a2.d[["model"]]$`as.factor(incumbent)`, test = "Satterthwaite")[,"SE"]
rank.num.a <- ame.p(a2.d, dv = "ranks_competed", iv = "ln_pq_land",seed = 1104, nsim = 1000)


##### -- Women -- #####
### --- PQs - Share --- ###
# remove office
summary(w2.c <- glm(ranks_competed ~ - 1+ as.factor(incumbent) + rank_before_prob + 
                      elected+ theta_mean_land_deviation + ln_pq_total+ pq_both_coverage_land ,  data = dat[dat$sex == "female",], family = poisson("log")))
w2.c.se.rob <- coeftest(w2.c, vcov = vcovHC(w2.c, "HC"))[,2];coeftest(w2.c, vcov = vcovHC(w2.c, "HC"))
w2.c.se.cl <- coef_test(w2.c, vcov = "CR2", 
                        cluster = w2.c[["model"]]$`as.factor(incumbent)`, test = "Satterthwaite")[,"SE"]
rank.cov.w <- ame.p(w2.c, dv = "ranks_competed", iv = "pq_both_coverage_land",seed = 1104, nsim = 1000)

### --- PQs - Number --- ###
summary(w2.d <- glm(ranks_competed ~ - 1+ as.factor(incumbent) + rank_before_prob + 
                      elected+ theta_mean_land_deviation + ln_pq_total + ln_pq_land,  data = dat[dat$sex == "female",], family = poisson("log")))
w2.d.se.rob <- coeftest(w2.d, vcov = vcovHC(w2.d, "HC"))[,2];coeftest(w2.d, vcov = vcovHC(w2.d, "HC"))
w2.d.se.cl <- coef_test(w2.d, vcov = "CR2", 
                        cluster = w2.d[["model"]]$`as.factor(incumbent)`, test = "Satterthwaite")[,"SE"]
rank.num.w <- ame.p(w2.d, dv = "ranks_competed", iv = "ln_pq_land",seed = 1104, nsim = 1000)



##### -- Men -- #####
### --- PQs - Share --- ###
# remove smd_nomination
summary(m2.c <- glm(ranks_competed ~ - 1+ as.factor(incumbent) + rank_before_prob +  
                      office + elected+ theta_mean_land_deviation+ ln_pq_total +pq_both_coverage_land ,  data = dat[dat$sex == "male",], family = poisson("log")))
m2.c.se.rob <- coeftest(m2.c, vcov = vcovHC(m2.c, "HC"))[,2];coeftest(m2.c, vcov = vcovHC(m2.c, "HC"))
m2.c.se.cl <- coef_test(m2.c, vcov = "CR2", 
                        cluster = m2.c[["model"]]$`as.factor(incumbent)`, test = "Satterthwaite")[,"SE"]
rank.cov.m <- ame.p(m2.c, dv = "ranks_competed", iv = "pq_both_coverage_land",seed = 1104, nsim = 1000)

### --- PQs - Number --- ###
summary(m2.d <- glm(ranks_competed ~ - 1+ as.factor(incumbent) + rank_before_prob + 
                      office + elected+ theta_mean_land_deviation +ln_pq_total + ln_pq_land,  data = dat[dat$sex == "male",], family = poisson("log")))
m2.d.se.rob <- coeftest(m2.d, vcov = vcovHC(m2.d, "HC"))[,2];coeftest(m2.d, vcov = vcovHC(m2.d, "HC"))
m2.d.se.cl <- coef_test(m2.d, vcov = "CR2", 
                        cluster = m2.d[["model"]]$`as.factor(incumbent)`, test = "Satterthwaite")[,"SE"]
rank.num.m <- ame.p(m2.d, dv = "ranks_competed", iv = "ln_pq_land",seed = 1104, nsim = 1000)


#################################################
# Table 7, Appendix
#################################################

stargazer(list( a2.c, a2.d, a2.c, a2.d),
          se = list( a2.c.se.rob, a2.d.se.rob, 
                     a2.c.se.cl, a2.d.se.cl), omit = c("incumbent"))

#################################################
# Table 10, Appendix
#################################################

stargazer(list( w2.c, m2.c, w2.d, m2.d,w2.c, m2.c, w2.d, m2.d),
          se = list( w2.c.se.rob, m2.c.se.rob, w2.d.se.rob, m2.d.se.rob,
                     w2.c.se.cl, m2.c.se.cl, w2.d.se.cl, m2.d.se.cl), omit = c("incumbent"))


dispersiontest(w2.c);dispersiontest(w2.d)
dispersiontest(m2.c);dispersiontest(m2.d)
#################################################
# 2) difference in rank
#################################################
dat$rank_diff <- dat$rank_before - dat$rank_final
dat$rank_diff[!is.na(dat$rank_diff) &  dat$rank_final == 0] <- -1

##### -- All -- #####
### --- PQs - Share --- ###
summary(a4.c <- lm(rank_diff ~ - 1+ as.factor(incumbent) + rank_before_prob + 
                     office + elected+ theta_mean_land_deviation+ ln_pq_total + pq_both_coverage_land ,  data = dat))
a4.c.se.rob <- coeftest(a4.c, vcov = vcovHC(a4.c, "HC"))[,2];coeftest(a4.c, vcov = vcovHC(a4.c, "HC"))
a4.c.se.cl <- coef_test(a4.c, vcov = "CR2", 
                        cluster = a4.c[["model"]]$`as.factor(incumbent)`, test = "Satterthwaite")[,"SE"]
diff.cov.a <- ame.p(a4.c, dv = "rank_diff", iv = "pq_both_coverage_land",seed = 1104, nsim = 1000)

### --- PQs - Number --- ###
summary(a4.d <- lm(rank_diff ~ - 1+ as.factor(incumbent) + rank_before_prob + 
                     office + elected+ theta_mean_land_deviation + ln_pq_total + ln_pq_land  ,  data = dat))
a4.d.se.rob <- coeftest(a4.d, vcov = vcovHC(a4.d, "HC"))[,2];coeftest(a4.d, vcov = vcovHC(a4.d, "HC"))
a4.d.se.cl <- coef_test(a4.d, vcov = "CR2", 
                        cluster = a4.d[["model"]]$`as.factor(incumbent)`, test = "Satterthwaite")[,"SE"]
diff.num.a <- ame.p(a4.d, dv = "rank_diff", iv = "ln_pq_land",seed = 1104, nsim = 1000)



##### -- Women -- #####
### --- PQs - Share --- ###
summary(w4.c <- lm(rank_diff ~ - 1+ as.factor(incumbent) + rank_before_prob + 
                     office + elected+ theta_mean_land_deviation+ ln_pq_total + pq_both_coverage_land ,  data = dat[dat$sex == "female",]))
w4.c.se.rob <- coeftest(w4.c, vcov = vcovHC(w4.c, "HC"))[,2];coeftest(w4.c, vcov = vcovHC(w4.c, "HC"))
w4.c.se.cl <- coef_test(w4.c, vcov = "CR2", 
                        cluster = w4.c[["model"]]$`as.factor(incumbent)`, test = "Satterthwaite")[,"SE"]
diff.cov.w <- ame.p(w4.c, dv = "rank_diff", iv = "pq_both_coverage_land",seed = 1104, nsim = 1000)

### --- PQs - Number --- ###
summary(w4.d <- lm(rank_diff ~ - 1+ as.factor(incumbent) + rank_before_prob + 
                     office + elected+ theta_mean_land_deviation + ln_pq_total + ln_pq_land  ,  data = dat[dat$sex == "female",]))
w4.d.se.rob <- coeftest(w4.d, vcov = vcovHC(w4.d, "HC"))[,2];coeftest(w4.d, vcov = vcovHC(w4.d, "HC"))
w4.d.se.cl <- coef_test(w4.d, vcov = "CR2", 
                        cluster = w4.d[["model"]]$`as.factor(incumbent)`, test = "Satterthwaite")[,"SE"]
diff.num.w <- ame.p(w4.d, dv = "rank_diff", iv = "ln_pq_land",seed = 1104, nsim = 1000)


##### -- Men -- #####
### --- PQs - Share --- ###
summary(m4.c <- lm(rank_diff ~ - 1+ as.factor(incumbent) + rank_before_prob + 
                     office + elected+ theta_mean_land_deviation + ln_pq_total + pq_both_coverage_land ,  data = dat[dat$sex == "male",]))#, family = binomial("logit")))
m4.c.se.rob <- coeftest(m4.c, vcov = vcovHC(m4.c, "HC"))[,2];coeftest(m4.c, vcov = vcovHC(m4.c, "HC"))
m4.c.se.cl <- coef_test(m4.c, vcov = "CR2", 
                        cluster = m4.c[["model"]]$`as.factor(incumbent)`, test = "Satterthwaite")[,"SE"]
diff.cov.m <- ame.p(m4.c, dv = "rank_diff", iv = "pq_both_coverage_land",seed = 1104, nsim = 1000)

### --- PQs - Number --- ###
summary(m4.d <- lm(rank_diff ~ - 1+ as.factor(incumbent) + rank_before_prob + 
                     office + elected+ theta_mean_land_deviation+ ln_pq_total + ln_pq_land  ,  data = dat[dat$sex == "male",]))#, family = binomial("logit")))
m4.d.se.rob <- coeftest(m4.d, vcov = vcovHC(m4.d, "HC"))[,2];coeftest(m4.d, vcov = vcovHC(m4.d, "HC"))
m4.d.se.cl <- coef_test(m4.d, vcov = "CR2", 
                        cluster = m4.d[["model"]]$`as.factor(incumbent)`, test = "Satterthwaite")[,"SE"]
diff.num.m <- ame.p(m4.d, dv = "rank_diff", iv = "ln_pq_land",seed = 1104, nsim = 1000)


#################################################
# Table 8, Appendix
#################################################

stargazer(list( a4.c, a4.d,  a4.c, a4.d),
          se = list( a4.c.se.rob, a4.d.se.rob, 
                     a4.c.se.cl, a4.d.se.cl), omit = c("incumbent"))

#################################################
# Table 11, Appendix
#################################################

stargazer(list( w4.c, m4.c, w4.d, m4.d,w4.c, m4.c, w4.d, m4.d),
          se = list( w4.c.se.rob, m4.c.se.rob, w4.d.se.rob, m4.d.se.rob,
                     w4.c.se.cl, m4.c.se.cl, w4.d.se.cl, m4.d.se.cl), omit = c("incumbent"))


#################################################
# 3) margin for final rank
#################################################

##### -- All -- #####
### --- PQs - Share --- ###
summary(a1.c <- lm(margin_final ~ -1 +  as.factor(incumbent) + rank_before_prob + 
                     office + elected+ theta_mean_land_deviation + ln_pq_total + pq_both_coverage_land ,  data =  dat))
a1.c.se.rob <- coeftest(a1.c, vcov = vcovHC(a1.c, "HC"))[,2];coeftest(a1.c, vcov = vcovHC(a1.c, "HC"))
a1.c.se.cl <- coef_test(a1.c, vcov = "CR2", 
                        cluster = a1.c[["model"]]$`as.factor(incumbent)`, test = "Satterthwaite")[,"SE"]
mar.cov.a <- ame.p(a1.c, dv = "margin_final", iv = "pq_both_coverage_land",seed = 1104, nsim = 1000)

### --- PQs - Number --- ###
summary(a1.d <- lm(margin_final ~ -1 +  as.factor(incumbent) + rank_before_prob +
                     office + elected+ theta_mean_land_deviation + ln_pq_total + ln_pq_land    ,  data = dat))
a1.d.se.rob <- coeftest(a1.d, vcov = vcovHC(a1.d, "HC"))[,2];coeftest(a1.d, vcov = vcovHC(a1.d, "HC"))
a1.d.se.cl <- coef_test(a1.d, vcov = "CR2", 
                        cluster = a1.d[["model"]]$`as.factor(incumbent)`, test = "Satterthwaite")[,"SE"]
mar.num.a <- ame.p(a1.d, dv = "margin_final", iv = "ln_pq_land",seed = 1104, nsim = 1000)


##### -- Women -- #####
### --- PQs - Share --- ###
summary(w1.c <- lm(margin_final ~ -1 +  as.factor(incumbent) + rank_before_prob + 
                     office + elected+ theta_mean_land_deviation + ln_pq_total + pq_both_coverage_land  ,  data =  dat[dat$sex == "female",]))
w1.c.se.rob <- coeftest(w1.c, vcov = vcovHC(w1.c, "HC"))[,2];coeftest(w1.c, vcov = vcovHC(w1.c, "HC"))
w1.c.se.cl <- coef_test(w1.c, vcov = "CR2", 
                        cluster = w1.c[["model"]]$`as.factor(incumbent)`, test = "Satterthwaite")[,"SE"]
mar.cov.w <- ame.p(w1.c, dv = "margin_final", iv = "pq_both_coverage_land",seed = 1104, nsim = 1000)

### --- PQs - Number --- ###
summary(w1.d <- lm(margin_final ~ -1 +  as.factor(incumbent) + rank_before_prob +
                     office + elected+ theta_mean_land_deviation + ln_pq_total + ln_pq_land   ,  data = dat[dat$sex == "female",]))
w1.d.se.rob <- coeftest(w1.d, vcov = vcovHC(w1.d, "HC"))[,2];coeftest(w1.d, vcov = vcovHC(w1.d, "HC"))
w1.d.se.cl <- coef_test(w1.d, vcov = "CR2", 
                        cluster = w1.d[["model"]]$`as.factor(incumbent)`, test = "Satterthwaite")[,"SE"]
mar.num.w <- ame.p(w1.d, dv = "margin_final", iv = "ln_pq_land",seed = 1104, nsim = 1000)


##### -- Men -- #####
### --- PQs - Share --- ###
summary(m1.c <- lm(margin_final ~ -1 +  as.factor(incumbent) + rank_before_prob + 
                     office + elected+ theta_mean_land_deviation + ln_pq_total+ pq_both_coverage_land   ,  data =  dat[dat$sex == "male",]))
m1.c.se.rob <- coeftest(m1.c, vcov = vcovHC(m1.c, "HC"))[,2];coeftest(m1.c, vcov = vcovHC(m1.c, "HC"))
m1.c.se.cl <- coef_test(m1.c, vcov = "CR2", 
                        cluster = m1.c[["model"]]$`as.factor(incumbent)`, test = "Satterthwaite")[,"SE"]
mar.cov.m <- ame.p(m1.c, dv = "margin_final", iv = "pq_both_coverage_land",seed = 1104, nsim = 1000)

### --- PQs - Number --- ###
summary(m1.d <- lm(margin_final ~ -1 +  as.factor(incumbent) + rank_before_prob +
                     office + elected+ ln_pq_total + ln_pq_land + theta_mean_land_deviation,  data = dat[dat$sex == "male",]))
m1.d.se.cl <- coef_test(m1.d, vcov = "CR2", 
                        cluster = m1.d[["model"]]$`as.factor(incumbent)`, test = "Satterthwaite")[,"SE"]
m1.d.se.rob <- coeftest(m1.d, vcov = vcovHC(m1.d, "HC"))[,2];coeftest(m1.d, vcov = vcovHC(m1.d, "HC"))

mar.num.m <- ame.p(m1.d, dv = "margin_final", iv = "ln_pq_land",seed = 1104, nsim = 1000)


#################################################
# Table 9, Appendix
#################################################

stargazer(list(  a1.c, a1.d,  a1.c, a1.d),
          se = list(  a1.c.se.rob, a1.d.se.rob, 
                      a1.c.se.cl, a1.d.se.cl), omit = c("incumbent"))

#################################################
# Table 12, Appendix
#################################################

stargazer(list(  w1.c, m1.c, w1.d, m1.d,w1.c, m1.c, w1.d, m1.d),
          se = list(  w1.c.se.rob, m1.c.se.rob, w1.d.se.rob, m1.d.se.rob,
                      w1.c.se.cl, m1.c.se.cl, w1.d.se.cl, m1.d.se.cl), omit = c("incumbent"))

#################################################
# 4) election probability of final rank
#################################################

##### -- All -- #####
### --- PQs - Share --- ###
summary(a1.c <- glm(rank_final_prob ~ -1 +as.factor(incumbent) + rank_before_prob + 
                      office + elected+ theta_mean_land_deviation + ln_pq_total + pq_both_coverage_land  , family = quasibinomial("logit"),  data =  dat))
a1.c.se.rob <- coeftest(a1.c, vcov = vcovHC(a1.c, "HC"))[,2];coeftest(a1.c, vcov = vcovHC(a1.c, "HC"))
a1.c.se.cl <- coef_test(a1.c, vcov = "CR2",  cluster = dat$land, test = "Satterthwaite")[,"SE"]

### --- PQs - Number --- ###
summary(a1.d <- glm(rank_final_prob ~ -1 +as.factor(incumbent) + rank_before_prob +
                      office + elected+ theta_mean_land_deviation + ln_pq_total + ln_pq_land   , family = quasibinomial("logit"), data = dat))
a1.d.se.rob <- coeftest(a1.d, vcov = vcovHC(a1.d, "HC"))[,2];coeftest(a1.d, vcov = vcovHC(a1.d, "HC"))
a1.d.se.cl <- coef_test(a1.d, vcov = "CR2",   cluster = dat$land, test = "Satterthwaite")[,"SE"]


##### -- Women -- #####
### --- PQs - Share --- ###
summary(w1.c <- glm(rank_final_prob ~ -1 +as.factor(incumbent) + rank_before_prob +
                      office + elected+ theta_mean_land_deviation + ln_pq_total + pq_both_coverage_land  , family = quasibinomial("logit"), data =  dat[dat$sex == "female",]))
w1.c.se.rob <- coeftest(w1.c, vcov = vcovHC(w1.c, "HC"))[,2];coeftest(w1.c, vcov = vcovHC(w1.c, "HC"))
w1.c.se.cl <- coef_test(w1.c, vcov = "CR2", cluster = dat$land[dat$sex == "female"], test = "Satterthwaite")[,"SE"]

### --- PQs - Number --- ###
summary(w1.d <- glm(rank_final_prob ~ -1 +as.factor(incumbent)+ rank_before_prob +
                      office + elected+ theta_mean_land_deviation + ln_pq_total + ln_pq_land   , family = quasibinomial("logit"),  data = dat[dat$sex == "female",]))
w1.d.se.rob <- coeftest(w1.d, vcov = vcovHC(w1.d, "HC"))[,2];coeftest(w1.d, vcov = vcovHC(w1.d, "HC"))
w1.d.se.cl <- coef_test(w1.d, vcov = "CR2",  cluster = dat$land[dat$sex == "female"], test = "Satterthwaite")[,"SE"]


##### -- Men -- #####
### --- PQs - Share --- ###
summary(m1.c <- glm(rank_final_prob ~ -1 +as.factor(incumbent) + rank_before_prob + 
                      office + elected+ theta_mean_land_deviation + ln_pq_total+ pq_both_coverage_land   , family = quasibinomial("logit"), data =  dat[dat$sex == "male",]))
m1.c.se.rob <- coeftest(m1.c, vcov = vcovHC(m1.c, "HC"))[,2];coeftest(m1.c, vcov = vcovHC(m1.c, "HC"))
m1.c.se.cl <- coef_test(m1.c, vcov = "CR2",  cluster = dat$land[dat$sex == "male"], test = "Satterthwaite")[,"SE"]


### --- PQs - Number --- ###
summary(m1.d <- glm(rank_final_prob ~ -1 +as.factor(incumbent) + rank_before_prob +
                      office + elected+ ln_pq_total + ln_pq_land + theta_mean_land_deviation, family = quasibinomial("logit"),  data = dat[dat$sex == "male",]))
m1.d.se.cl <- coef_test(m1.d, vcov = "CR2",  cluster = dat$land[dat$sex == "male"], test = "Satterthwaite")[,"SE"]
m1.d.se.rob <- coeftest(m1.d, vcov = vcovHC(m1.d, "HC"))[,2];coeftest(m1.d, vcov = vcovHC(m1.d, "HC"))


#################################################
# Table 14, Appendix
#################################################

stargazer(list(  w1.c, m1.c, w1.d, m1.d,w1.c, m1.c, w1.d, m1.d),
          se = list(  w1.c.se.rob, m1.c.se.rob, w1.d.se.rob, m1.d.se.rob,
                      w1.c.se.cl, m1.c.se.cl, w1.d.se.cl, m1.d.se.cl), omit = "incumbent")

 
#################################################
# Table 5, Appendix
#################################################

stargazer(dat[!is.na(dat$theta_mean_land_deviation) & !is.na(dat$rank_diff),
              c("rank_diff","margin_final", "ranks_competed",  
                "pq_both_coverage_land","pq_both_num_land" , "pq_total",  
                "rank_before_prob","mandates_last","government","age",
                "smd_nomination","office","elected","theta_mean_land_deviation")])


#################################################
# Figure 4, Manuscript
#################################################

pData4 <- data.frame(rbind(
  cbind(outcome = "No. of Ranks Competed"  , iv = "Share of \nKreise Covered", ame = rank.cov.a$ame95[1], ci = 95, lower = rank.cov.a$ame95[2], upper = rank.cov.a$ame95[3] ),
  cbind(outcome = "No. of Ranks Competed"  , iv = "Share of \nKreise Covered", ame = rank.cov.a$ame90[1], ci = 90, lower = rank.cov.a$ame90[2], upper = rank.cov.a$ame90[3] ),
  cbind(outcome = "No. of Ranks Competed"  , iv = "Number of \nLand Locations", ame = rank.num.a$ame95[1], ci = 95, lower = rank.num.a$ame95[2], upper = rank.num.a$ame95[3] ),
  cbind(outcome = "No. of Ranks Competed"  , iv = "Number of \nLand Locations", ame = rank.num.a$ame90[1], ci = 90, lower = rank.num.a$ame90[2], upper = rank.num.a$ame90[3] ),
  cbind(outcome = "Difference in Rank Position"  , iv = "Share of \nKreise Covered", ame = diff.cov.a$ame95[1], ci = 95, lower = diff.cov.a$ame95[2], upper = diff.cov.a$ame95[3] ),
  cbind(outcome = "Difference in Rank Position"  , iv = "Share of \nKreise Covered", ame = diff.cov.a$ame90[1], ci = 90, lower = diff.cov.a$ame90[2], upper = diff.cov.a$ame90[3] ),
  cbind(outcome = "Difference in Rank Position"  , iv = "Number of \nLand Locations", ame = diff.num.a$ame95[1], ci = 95, lower = diff.num.a$ame95[2], upper = diff.num.a$ame95[3] ),
  cbind(outcome = "Difference in Rank Position"  , iv = "Number of \nLand Locations", ame = diff.num.a$ame90[1], ci = 90, lower = diff.num.a$ame90[2], upper = diff.num.a$ame90[3] ),
  cbind(outcome = "Vote Margin for Final Rank"  , iv = "Share of \nKreise Covered", ame = mar.cov.a$ame95[1], ci = 95, lower = mar.cov.a$ame95[2], upper = mar.cov.a$ame95[3] ),
  cbind(outcome = "Vote Margin for Final Rank"  , iv = "Share of \nKreise Covered", ame = mar.cov.a$ame90[1], ci = 90, lower = mar.cov.a$ame90[2], upper = mar.cov.a$ame90[3] ),
  cbind(outcome = "Vote Margin for Final Rank"  , iv = "Number of \nLand Locations", ame = mar.num.a$ame95[1], ci = 95, lower = mar.num.a$ame95[2], upper = mar.num.a$ame95[3] ),
  cbind(outcome = "Vote Margin for Final Rank"  , iv = "Number of \nLand Locations", ame = mar.num.a$ame90[1], ci = 90, lower = mar.num.a$ame90[2], upper = mar.num.a$ame90[3] )))


row.names(pData4) <- 1:nrow(pData4)
pData4$ame <- as.numeric(as.character(pData4$ame))
pData4$lower <- as.numeric(as.character(pData4$lower))
pData4$upper <- as.numeric(as.character(pData4$upper))
pData4$ci <- as.numeric(as.character(pData4$ci))
pData4$outcome <- as.character(pData4$outcome)
pData4$outcome_f <- factor(pData4$outcome, levels = c("No. of Ranks Competed","Vote Margin for Final Rank","Difference in Rank Position"))


library(ggplot2)

p = ggplot(pData4[pData4$ci == 95 & pData4$iv != "Deviation from\nLandesverband",])
p = p + geom_pointrange(data = pData4[pData4$ci == 90 & pData4$iv != "Deviation from\nLandesverband",], aes(x = iv, y = ame, ymin = lower, ymax = upper))
p = p + geom_pointrange(data = pData4[pData4$ci == 95 & pData4$iv != "Deviation from\nLandesverband",],aes(x = iv, y = ame, ymin = lower, ymax = upper), lty = 2)
p = p + geom_hline(aes(yintercept = 0), linetype = 2, lwd = 1.05) + xlab("") + ylab("AME")
p = p + facet_wrap(vars(outcome_f), scales = 'free_x')
p = p + coord_flip() + theme(axis.text = element_text(size = rel(1), face = "bold", colour = "black"),
                             axis.ticks = element_line(colour = "black"),
                             legend.key = element_rect(fill = "white"),
                             plot.margin = unit(c(0, 0, 0, 0), "cm"),
                             panel.background = element_rect(fill = "white", colour = NA),
                             panel.border = element_rect(fill = NA, colour = "grey50"),
                             strip.background = element_rect(fill = "white", colour = "grey50"),
                             strip.text.x = element_text(size = 12, colour = "black", face = "bold"),
                             strip.text.y = element_text(size = 10, colour = "black", face = "bold"))  
p


#################################################
# Figure 5, Manuscript
#################################################

pData5 <- data.frame(rbind(
  cbind(gender = "Woman",outcome = "No. of Ranks Competed"  , iv = "Share of \nKreise Covered", ame = rank.cov.w$ame95[1], ci = 95, lower = rank.cov.w$ame95[2], upper = rank.cov.w$ame95[3] ),
  cbind(gender = "Woman",outcome = "No. of Ranks Competed"  , iv = "Share of \nKreise Covered", ame = rank.cov.w$ame90[1], ci = 90, lower = rank.cov.w$ame90[2], upper = rank.cov.w$ame90[3] ),
  cbind(gender = "Woman",outcome = "No. of Ranks Competed"  , iv = "Number of \nLand Locations", ame = rank.num.w$ame95[1], ci = 95, lower = rank.num.w$ame95[2], upper = rank.num.w$ame95[3] ),
  cbind(gender = "Woman",outcome = "No. of Ranks Competed"  , iv = "Number of \nLand Locations", ame = rank.num.w$ame90[1], ci = 90, lower = rank.num.w$ame90[2], upper = rank.num.w$ame90[3] ),
  cbind(gender = "Woman",outcome = "Difference in Rank Position"  , iv = "Share of \nKreise Covered", ame = diff.cov.w$ame95[1], ci = 95, lower = diff.cov.w$ame95[2], upper = diff.cov.w$ame95[3] ),
  cbind(gender = "Woman",outcome = "Difference in Rank Position"  , iv = "Share of \nKreise Covered", ame = diff.cov.w$ame90[1], ci = 90, lower = diff.cov.w$ame90[2], upper = diff.cov.w$ame90[3] ),
  cbind(gender = "Woman",outcome = "Difference in Rank Position"  , iv = "Number of \nLand Locations", ame = diff.num.w$ame95[1], ci = 95, lower = diff.num.w$ame95[2], upper = diff.num.w$ame95[3] ),
  cbind(gender = "Woman",outcome = "Difference in Rank Position"  , iv = "Number of \nLand Locations", ame = diff.num.w$ame90[1], ci = 90, lower = diff.num.w$ame90[2], upper = diff.num.w$ame90[3] ),
  cbind(gender = "Woman",outcome = "Vote Margin for Final Rank"  , iv = "Share of \nKreise Covered", ame = mar.cov.w$ame95[1], ci = 95, lower = mar.cov.w$ame95[2], upper = mar.cov.w$ame95[3] ),
  cbind(gender = "Woman",outcome = "Vote Margin for Final Rank"  , iv = "Share of \nKreise Covered", ame = mar.cov.w$ame90[1], ci = 90, lower = mar.cov.w$ame90[2], upper = mar.cov.w$ame90[3] ),
  cbind(gender = "Woman",outcome = "Vote Margin for Final Rank"  , iv = "Number of \nLand Locations", ame = mar.num.w$ame95[1], ci = 95, lower = mar.num.w$ame95[2], upper = mar.num.w$ame95[3] ),
  cbind(gender = "Woman",outcome = "Vote Margin for Final Rank"  , iv = "Number of \nLand Locations", ame = mar.num.w$ame90[1], ci = 90, lower = mar.num.w$ame90[2], upper = mar.num.w$ame90[3] ),
  cbind(gender = "Man", outcome = "No. of Ranks Competed"  , iv = "Share of \nKreise Covered", ame = rank.cov.m$ame95[1], ci = 95, lower = rank.cov.m$ame95[2], upper = rank.cov.m$ame95[3] ),
  cbind(gender = "Man", outcome = "No. of Ranks Competed"  , iv = "Share of \nKreise Covered", ame = rank.cov.m$ame90[1], ci = 90, lower = rank.cov.m$ame90[2], upper = rank.cov.m$ame90[3] ),
  cbind(gender = "Man", outcome = "No. of Ranks Competed"  , iv = "Number of \nLand Locations", ame = rank.num.m$ame95[1], ci = 95, lower = rank.num.m$ame95[2], upper = rank.num.m$ame95[3] ),
  cbind(gender = "Man", outcome = "No. of Ranks Competed"  , iv = "Number of \nLand Locations", ame = rank.num.m$ame90[1], ci = 90, lower = rank.num.m$ame90[2], upper = rank.num.m$ame90[3] ),
  cbind(gender = "Man", outcome = "Difference in Rank Position"  , iv = "Share of \nKreise Covered", ame = diff.cov.m$ame95[1], ci = 95, lower = diff.cov.m$ame95[2], upper = diff.cov.m$ame95[3] ),
  cbind(gender = "Man", outcome = "Difference in Rank Position"  , iv = "Share of \nKreise Covered", ame = diff.cov.m$ame90[1], ci = 90, lower = diff.cov.m$ame90[2], upper = diff.cov.m$ame90[3] ),
  cbind(gender = "Man", outcome = "Difference in Rank Position"  , iv = "Number of \nLand Locations", ame = diff.num.m$ame95[1], ci = 95, lower = diff.num.m$ame95[2], upper = diff.num.m$ame95[3] ),
  cbind(gender = "Man", outcome = "Difference in Rank Position"  , iv = "Number of \nLand Locations", ame = diff.num.m$ame90[1], ci = 90, lower = diff.num.m$ame90[2], upper = diff.num.m$ame90[3] ),
  cbind(gender = "Man", outcome = "Vote Margin for Final Rank"  , iv = "Share of \nKreise Covered", ame = mar.cov.m$ame95[1], ci = 95, lower = mar.cov.m$ame95[2], upper = mar.cov.m$ame95[3] ),
  cbind(gender = "Man", outcome = "Vote Margin for Final Rank"  , iv = "Share of \nKreise Covered", ame = mar.cov.m$ame90[1], ci = 90, lower = mar.cov.m$ame90[2], upper = mar.cov.m$ame90[3] ),
  cbind(gender = "Man", outcome = "Vote Margin for Final Rank"  , iv = "Number of \nLand Locations", ame = mar.num.m$ame95[1], ci = 95, lower = mar.num.m$ame95[2], upper = mar.num.m$ame95[3] ),
  cbind(gender = "Man", outcome = "Vote Margin for Final Rank"  , iv = "Number of \nLand Locations", ame = mar.num.m$ame90[1], ci = 90, lower = mar.num.m$ame90[2], upper = mar.num.m$ame90[3] )))

row.names(pData5) <- 1:nrow(pData5)
pData5$ame <- as.numeric(as.character(pData5$ame))
pData5$lower <- as.numeric(as.character(pData5$lower))
pData5$upper <- as.numeric(as.character(pData5$upper))
pData5$ci <- as.numeric(as.character(pData5$ci))
pData5$outcome <- as.character(pData5$outcome)
pData5$outcome_f <- factor(pData5$outcome, levels = c("No. of Ranks Competed","Vote Margin for Final Rank","Difference in Rank Position")) #, "Election Probability of Final Rank"))


library(ggplot2)

p = ggplot(pData5[pData5$ci == 95 & pData5$iv != "Deviation from\nLandesverband" ,])
p = p + geom_pointrange(data = pData5[pData5$ci == 90& pData5$iv != "Deviation from\nLandesverband",], 
                        aes(x = iv, y = ame, ymin = lower, ymax = upper),
                        position=position_nudge(x = -0.05), shape = 16)
p = p + geom_pointrange(data = pData5[pData5$ci == 95 & pData5$iv != "Deviation from\nLandesverband",], 
                        aes(x = iv, y = ame, ymin = lower, ymax = upper),   lty = 2,
                        position=position_nudge(x = -0.05), shape = 16)
p = p + geom_hline(aes(yintercept = 0), linetype = 2) + xlab("") + ylab("AME")
p = p + facet_grid(gender ~ outcome_f, scales = 'free')
p = p + coord_flip() + theme(axis.text = element_text(size = rel(1), face = "bold", colour = "black"),
                             axis.ticks = element_line(colour = "black"),
                             legend.key = element_rect(fill = "white"),
                             plot.margin = unit(c(0, 0, 0, 0), "cm"),
                             panel.background = element_rect(fill = "white", colour = NA),
                             panel.border = element_rect(fill = NA, colour = "grey50"),
                             strip.background = element_rect(fill = "white", colour = "grey50"),
                             strip.text.x = element_text(size = 12, colour = "black", face = "bold"),
                             strip.text.y = element_text(size = 10, colour = "black", face = "bold"))  
p



#################################################
# explaining local responsiveness
#################################################
 
summary(loc.1 <- lm(pq_both_coverage_land ~ as.factor(land)+ government + rank_before_prob +  age +  office + election + elected + ln_pq_total + margin_final_before + sex
                    , data = dat))
loc.1.se.rob <- coeftest(loc.1, vcov = vcovHC(loc.1, "HC"))[,2];coeftest(loc.1, vcov = vcovHC(loc.1, "HC"))
loc.1.se.cl <- coef_test(loc.1, vcov = "CR2", 
                         cluster = loc.1[["model"]]$`as.factor(land)`, test = "Satterthwaite")[,"SE"]


summary(loc.2 <- lm(ln_pq_land ~ as.factor(land)+ government + rank_before_prob + age +  office +  election + elected + ln_pq_total + margin_final_before + sex 
                    , data = dat))
loc.2.se.rob <- coeftest(loc.2, vcov = vcovHC(loc.2, "HC"))[,2];coeftest(loc.2, vcov = vcovHC(loc.2, "HC"))
loc.2.se.cl <- coef_test(loc.2, vcov = "CR2", 
                         cluster = loc.2[["model"]]$`as.factor(land)`, test = "Satterthwaite")[,"SE"]


summary(loc.3 <- lm(pq_both_coverage_land ~ as.factor(land)+ government + rank_before_prob +  age +  office + election + elected + ln_pq_total + enec_final_before + sex   
                    , data = dat))
loc.3.se.rob <- coeftest(loc.3, vcov = vcovHC(loc.3, "HC"))[,2];coeftest(loc.3, vcov = vcovHC(loc.3, "HC"))
loc.3.se.cl <- coef_test(loc.3, vcov = "CR2", 
                         cluster = loc.3[["model"]]$`as.factor(land)`, test = "Satterthwaite")[,"SE"]


summary(loc.4 <- lm(ln_pq_land ~ as.factor(land)+ government + rank_before_prob + age +  office +  election + elected + ln_pq_total + enec_final_before + sex 
                    , data = dat))
loc.4.se.rob <- coeftest(loc.4, vcov = vcovHC(loc.4, "HC"))[,2];coeftest(loc.4, vcov = vcovHC(loc.4, "HC"))
loc.4.se.cl <- coef_test(loc.4, vcov = "CR2", 
                         cluster = loc.4[["model"]]$`as.factor(land)`, test = "Satterthwaite")[,"SE"]


summary(loc.5 <- lm(pq_both_coverage_land ~ as.factor(land)+ government + rank_before_prob +  age +  office + election + elected + ln_pq_total + num_cand_before1st + sex   
                    , data = dat))
loc.5.se.rob <- coeftest(loc.5, vcov = vcovHC(loc.5, "HC"))[,2];coeftest(loc.5, vcov = vcovHC(loc.5, "HC"))
loc.5.se.cl <- coef_test(loc.5, vcov = "CR2", 
                         cluster = loc.5[["model"]]$`as.factor(land)`, test = "Satterthwaite")[,"SE"]


summary(loc.6 <- lm(ln_pq_land ~ as.factor(land)+ government + rank_before_prob + age +  office +  election + elected + ln_pq_total + num_cand_before1st + sex 
                    , data = dat))
loc.6.se.rob <- coeftest(loc.6, vcov = vcovHC(loc.6, "HC"))[,2];coeftest(loc.6, vcov = vcovHC(loc.6, "HC"))
loc.6.se.cl <- coef_test(loc.6, vcov = "CR2", 
                         cluster = loc.6[["model"]]$`as.factor(land)`, test = "Satterthwaite")[,"SE"]


#################################################
# Table 12, Appendix
#################################################

stargazer(list(loc.1, loc.1, loc.2, loc.2,loc.3, loc.3, loc.4, loc.4,loc.5, loc.5, loc.6, loc.6),
          se = list(loc.1.se.rob, loc.1.se.cl, loc.2.se.rob, loc.2.se.cl,
                    loc.3.se.rob, loc.3.se.cl, loc.4.se.rob, loc.4.se.cl,
                    loc.5.se.rob, loc.5.se.cl, loc.6.se.rob, loc.6.se.cl),
          omit = c("land")
          )



stargazer(list(loc.1, loc.1, loc.2, loc.2,loc.3, loc.3),
          se = list(loc.1.se.rob, loc.1.se.cl, loc.2.se.rob, loc.2.se.cl,
                    loc.3.se.rob, loc.3.se.cl),
          omit = c("land"), out = "table12_1.rtf", type= "html"
)

stargazer(list(loc.4, loc.4,loc.5, loc.5, loc.6, loc.6),
          se = list( loc.4.se.rob, loc.4.se.cl,
                    loc.5.se.rob, loc.5.se.cl, loc.6.se.rob, loc.6.se.cl),
          omit = c("land"), out = "table12_2.rtf", type= "html"
)

#################################################
# Table 13, Appendix
#################################################

stargazer(list(loc.1, loc.1, loc.2, loc.2,loc.3, loc.3, loc.4, loc.4,loc.5, loc.5, loc.6, loc.6),
          se = list(loc.1.se.rob, loc.1.se.cl, loc.2.se.rob, loc.2.se.cl,
                    loc.3.se.rob, loc.3.se.cl, loc.4.se.rob, loc.4.se.cl,
                    loc.5.se.rob, loc.5.se.cl, loc.6.se.rob, loc.6.se.cl),
          omit = "land")



#################################################
###### Placebo Models ###### 
#################################################

##### -- Women -- #####
### --- PQs - Number --- ###
summary(w2.d <- glm(ranks_competed ~ -1+as.factor(incumbent)+  rank_before_prob + 
                      elected+ theta_mean_land_deviation + ln_pq_total  + ln_pq_other,  data = dat[dat$sex == "female",], family = poisson("log")))
w2.d.se.rob <- coeftest(w2.d, vcov = vcovHC(w2.d, "HC"))[,2];coeftest(w2.d, vcov = vcovHC(w2.d, "HC"))
w2.d.se.cl <- coef_test(w2.d, vcov = "CR2",  cluster = dat$incumbent[dat$sex == "female"], test = "Satterthwaite")[,"SE"]

##### -- Men -- #####
### --- PQs - Number --- ###
summary(m2.d <- glm(ranks_competed ~ -1+as.factor(incumbent) + rank_before_prob + 
                      office + elected+ theta_mean_land_deviation +ln_pq_total +  ln_pq_other,  data = dat[dat$sex == "male",], family = poisson("log")))
m2.d.se.rob <- coeftest(m2.d, vcov = vcovHC(m2.d, "HC"))[,2];coeftest(m2.d, vcov = vcovHC(m2.d, "HC"))
m2.d.se.cl <- coef_test(m2.d, vcov = "CR2",  cluster = dat$incumbent[dat$sex == "male"], test = "Satterthwaite")[,"SE"]

#################################################
# 2) difference in rank
#################################################

##### -- Women -- #####
### --- PQs - Number --- ###
summary(w4.d <- lm(rank_diff ~ -1+as.factor(incumbent) + rank_before_prob + 
                     office + elected+ theta_mean_land_deviation + ln_pq_total +  ln_pq_other,  data = dat[dat$sex == "female",]))
w4.d.se.rob <- coeftest(w4.d, vcov = vcovHC(w4.d, "HC"))[,2];coeftest(w4.d, vcov = vcovHC(w4.d, "HC"))
w4.d.se.cl <- coef_test(w4.d, vcov = "CR2", 
                        cluster = dat$incumbent[dat$sex == "female"], test = "Satterthwaite")[,"SE"]

##### -- Men -- #####
### --- PQs - Number --- ###
summary(m4.d <- lm(rank_diff ~ -1+as.factor(incumbent) + rank_before_prob + 
                     office + elected+ theta_mean_land_deviation+ ln_pq_total +  ln_pq_other,  data = dat[dat$sex == "male",]))#, family = binomial("logit")))
m4.d.se.rob <- coeftest(m4.d, vcov = vcovHC(m4.d, "HC"))[,2];coeftest(m4.d, vcov = vcovHC(m4.d, "HC"))
m4.d.se.cl <- coef_test(m4.d, vcov = "CR2",    cluster = dat$incumbent[dat$sex == "male"], test = "Satterthwaite")[,"SE"]


#################################################
# 3) margin for final rank
#################################################

##### -- Women -- #####
### --- PQs - Number --- ###
summary(w1.d <- lm(margin_final ~ -1+as.factor(incumbent) + rank_before_prob +
                     office + elected+ theta_mean_land_deviation + ln_pq_total +  ln_pq_other,  data = dat[dat$sex == "female",]))
w1.d.se.rob <- coeftest(w1.d, vcov = vcovHC(w1.d, "HC"))[,2];coeftest(w1.d, vcov = vcovHC(w1.d, "HC"))
w1.d.se.cl <- coef_test(w1.d, vcov = "CR2",    cluster = dat$incumbent[dat$sex == "female"], test = "Satterthwaite")[,"SE"]


##### -- Men -- #####
### --- PQs - Number --- ###
summary(m1.d <- lm(margin_final ~ -1+as.factor(incumbent) + rank_before_prob +
                     office + elected+ theta_mean_land_deviation + ln_pq_total +  ln_pq_other,  data = dat[dat$sex == "male",]))
m1.d.se.cl <- coef_test(m1.d, vcov = "CR2",  cluster = dat$incumbent[dat$sex == "male"], test = "Satterthwaite")[,"SE"]
m1.d.se.rob <- coeftest(m1.d, vcov = vcovHC(m1.d, "HC"))[,2];coeftest(m1.d, vcov = vcovHC(m1.d, "HC"))


#################################################
# Table 15, Appendix
#################################################

stargazer(list(  w2.d, m2.d, w4.d, m4.d,w1.d, m1.d,
                 w2.d, m2.d, w4.d, m4.d,w1.d, m1.d),
          se = list( w2.d.se.rob, m2.d.se.rob, w4.d.se.rob, m4.d.se.rob,w1.d.se.rob, m1.d.se.rob,
                     w2.d.se.cl, m2.d.se.cl, w4.d.se.cl, m4.d.se.cl,w1.d.se.cl, m1.d.se.cl), omit = "incumbent")



#################################################
# Figure 12, Appendix: Distribution of PQs
#################################################

par(mfrow = c(2,1),
    oma = c(3,5,1,0) + 0.1,
    mar = c(3,0,1,1) + 0.1)
hist(dat$pq_total, breaks = 50,  xlab = "", main = "Total No. of PQs", xaxt = "n", yaxt = "n", col = "white")
axis(1, seq(0,750,100), seq(0,750,100), col.ticks = "black", col = "white")
axis(2, seq(0,180,20), seq(0,180,20), col.ticks = "black", col = "white")

hist(dat$pq_both_num_land, breaks = 50,  xlab = "", main = "No. of Land Locations Mentioned in PQs", xaxt = "n", yaxt = "n", col = "white")
axis(1, seq(0,250,50), seq(0,250,50), col.ticks = "black", col = "white")
axis(2, seq(0,200,20), seq(0,200,20), col.ticks = "black", col = "white")

